Migration of a subsurface wavefield in reflection seismics: A mathematical study 
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ABSTRACT 

In this pedagogically motivated work, the process of migration in reflection seismics has been consid- 
ered from a rigorously mathematical viewpoint. An inclined subsurface reflector with a constant dipping 
angle has been shown to cause a shift in the normal moveout equation, with the peak of the moveout curve 
tracing an elliptic locus. Since any subsurface reflector actually has a non-uniform spatial variation, the 
use of a more comprehensive principle of migration, by adopting the wave equation, has been argued to be 
necessary. By this approach an expression has been derived for both the amplitude and the phase of a sub- 
surface wavefield with vertical velocity variation. This treatment has entailed the application of the WKB 
approximation, whose self-consistency has been established by the fact that the logarithmic variation of 
the velocity is very slow in the vertical direction, a feature that is much more strongly upheld at increas- 
ingly greater subsurface depths. Finally, it has been demonstrated that for a planar subsurface wavefield, 
there is an equivalence between the constant velocity Stolt Migration algorithm and the stationary phase 
approximation method (by which the origin of the reflected subsurface signals is determined). 

Subject headings: reflection seismics — migration — mathematical methods 



1. Reflection seismology: A basic introduction 

In seismic processing, migration of seismic data is an exercise of paramount importance. Unless this aspect of 
processing is addressed satisfactorily, the possibility of making misleading estimates of subsurface depths remains 
wide open. Even by a simple quantitative estimate it can be shown that i n many cases the error in the spatial mea- 



surement of the true reflector position could be of the order of a kilometre ( YilmazfeoOl ). Naturally, when it comes to 
exploring hydrocarbon reserves, this lack of precision will translate into a huge and wasted monetary expense. So an 
accurate representation of the depth and the local dip of subsurfa ce reflectors is impe rative. Migration of seismic data 
goes a long way in resolving this most difficult question ( Lowriefl997 : YilmazlEoOl ). 



It should be useful f or any study on the relevance o f migration to start with a simple and pedagogic understanding 
of reflection seismology dLowrielll997 ; Mari et aO 19991) . whose primary purpose is to find the depths to reflecting in- 



terfaces and the corresponding seismic velocities of subsurface rock layers. This process first necessitates the artificial 
generation of mechanical waves at predetermined times and spatial positions on the surface of the earth. These waves 
pass through various earth layers and are reflected from the subsurface interfaces back to the surface of the earth, 
where they are recorded by receivers. For immediate purposes, the waves of intere st are compressional (longitudinal) 
waves, which are also referred to as primary waves ( Lowriefl997 ; Mari et aO 19991) . 



Considering the simplest possible case of this system, one might suppose the surface of the earth to be absolutely 
flat. From this surface a single source is creating waves which travel through the subsurface and are subsequently 
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Fig. 1. — A schematic representation of reflection seismics from a flat subsurface stratum, using a single source-receiver pair. 
The waves originate from the source position, S. They are reflected from an absolutely flat interface in the subsurface, back to the 
receiver position, R. The midpoint between S and R, shown as M, is vertically above the point from where the reflection occurs, 
indicated as D. For many source-receiver arrays placed as symmetrically about M, as it has been shown in the figure, the reflection 
point D would be a common depth point for all such source-receiver pairs. The source-to-receiver distance, SR, is x. The depth to 
the reflecting interface, MD, is given as d. The signal takes a time t, travelling with a constant velocity v, to reach R from S, after 
undergoing a reflection at D. Triangle SDR is an isosceles triangle, while triangle MDR is a right-angled triangle. 



reflected back to the surface from a single stratum that is also equally flat. This has been schematically represented in 
Fig.Q] In this simple situation it is possible to set down the geometrical locus for the physical raypath in the x—t plane, 
with x defining the source-to-receiver position, and with the receiver itself gathering reflected signals after an interval 
of time, t. 



To that end the Pythagoras theorem may be applied to the triangle MDR in Fig. [T] under the prescription that 
MR = x/2, MD = d (the vertical depth of the reflecting interface from the surface of the earth), and DR = vt/2 (with 
v being the constant velocity of the waves propagating through the subsurface). Followed by some simple algebraic 
manipulations, this exercise will deliver the x-t locus corresponding to the actual raypath, known also as the normal 
moveout equation ( Lowriefl997 ). as 

t J \VtoJ 



(1) 



which is an expression that very recognisably bears the canonical form of a hyperbola, with the parameter to having 
been defined as to = 2d/v. The solution given by Eq. (Q]) has been geometrically depicted by the continuous curve in 
Fig. |2] From the plot it may also be appreciated that different values of d, as a parameter in Eq. (Q3, will yield a family 
of hyperbolae (assuming that v retains the same constant value all through). The turning points, the maxima in this 
case, of this family of hyperbolae (represented by the single continuous curve in Fig. [2]), obtained under the condition 
dt/dx = 0, will correspond to the coordinate (0, to) m the x-t plane. And for any fixed value of x in the x-t plane, the 
entire range of points along the t axis (with t customarily having been scaled along the vertical direction), c utting down 
the en tire family of hyperbolae, will define what is known in seismic processing terms as a seismic trace jMari et al.l 



1999b . 

Recast in a slightly different form, Eq. (HJ will read as 



\ 1+|p> ( 2 ) 
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Fig. 2. — The locus of the raypath between a source and a receiver in the x—t plane. In this plot, t increases in the downward 
direction along the vertical axis. The origin of coordinates, (0, 0), is at O. The continuous curve corresponds to the source-receiver 
array shown in Fig.Q] The peak of this curve is at the point P, whose coordinates are (0, to) in this x-t plot. The dotted curve 
corresponds to the situation shown in Fig. [5] The two horizontal lines are the tangents passing through the maxima of the two 
hyperbolae. It is evident from this plot that inclined subsurface reflectors cause a shift in the locus of the raypath. In this case the 
maximum point on the curve has shifted to the point Pi. With continuous change in the dip angle of the subsurface reflector, the 
path traced out by the peaks of the hyperbolae is an ellipse. 



in which to can physically be seen to represent the "echo time" — the two-way travel time for the wave, corresponding 
to a vertical reflection from the subsurface reflector. The term under the square root sign is the normal moveout factor, 
which appears because the wave reachi ng the receive r at a distance x from the shot point (i.e. the location of the 
source), has not been reflected vertically (Lowri el 1997b . Under the condition that the receiver distance, x, is much less 
than the reflector depth, d, i.e. x <C d, it is a matter of some simple algebraic steps to arrive at an approximate solution 
(through a binomial expansion) for the normal moveout, At = t — to, as 

2 

2v z to 

The primary objective here is to determine the depth, d, of the reflector. The vertical echo time, to, and the normal 
moveout time, At, can be obtained from the reflection data. The corresponding receiver distance, x, will also be known 
from the data. Taken together, all of these will allow for determining the value of v from Eq. (0. This, with the help 
of the definition of to, will subsequently also allow for the determination of d. In all of this there is a clear signal 
that a precise extraction of the value of v from the reflection data is a principal necessity in determining the reflector 
depth accurately. Even the simplest possible case that is being studied here emphasises this fact. Knowing the precise 
value of the velocity will be an equally important issue when one considers more complicated instances of reflection 
seismology involving non-horizontal reflectors of continually varying grad ients. Besides this, with v being dependen t 
on the elastic properties of the material through which the wave propagates ( Landau & Lifshitzll986 ; Mari et aljl999 ), 



getting a correct measure of the velocity also conveys an impression of the true nature of the subsurface material. 



2. Dipping reflectors: A case for migration 



It is scarcely to be expected that the layering of actual geological strata will conform to the neat and orderly 
picture implied by Fig. Q] If anything, not only will the subsurface strata be quite great in number (instead of the 
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Fig. 3. — A schematic representation of reflection seismics from an inclined subsurface reflector, using a single source-receiver 
pair placed symmetrically about the midpoint, M. The inclined reflector is at a dip angle 9 with respect to the horizontal. The 
actual point of subsurface reflection, Di, is not vertically below M (unlike the point D), and its location will vary for each set of 
source-receiver pair placed symmetrically about M. There will be no common depth point for these multiple source-receiver pairs, 
because of the asymmetry arising from the dip. The perpendicular distance from the source to the reflector, SN = d. The total 
travel path for the wave is SDi + DiR = vt, and the source-to-receiver distance, SR = x. 



single stratum shown in Fig. [TJ, but also each stratum will have complex gradients in space. This will cause for much 
difficulty in devising an accurate theoretical model for the physics of reflection seismics. This fact can be amply 
illustrated even by considering the very simple case of a single dipping interface with a constant spatial gradient. 

When the reflecting interface is inclined at a constant angle 9, with respect to the horizontal, as it has been 
shown in Fig. [3] it will be necessary to make suitable modifications in the analysis presented in Section Q] to retain 
the canonical form of the hyperbolic normal moveout equation that Eq. (Q]i represents. This can indeed be achieved 
by noting in Fig. [3] that the perpendicular distance from the source to the reflector, SN = d, the source-to-receiver 
distance, SR = x, and the total travel path, SDi + DiR = vt. With this information, even upon accounting for the 
dip in the reflector, the invariant canonical form of the hyperbola can be obtained as 




with x and to being defined by the transformations, x = x — 2d sin 9 and to = (2d cos 9) / v, respectively. 

From Eq. (0J, the coordinates of the maximum point of the hyperbola (x m , t m ) in the x—t plane (with t increasing 
in the downward direction) will be given by the condition dt/dx = 0. This will place the peak of the hyperbola at 
(— 2dsin9, to). Going back to the case of the horizontal reflector (9 = 0), given by Eq. (Q}, the peak of the hyperbola 
here is to be seen at (0, to). Therefore, a distinct shift is seen to arise because of the dip of the subsurface reflector, 
and this shift has been shown by the position of the dotted hyperbola in Fig. [2] 

The locus of this shift can also be traced on the x—t plane. Noting that the coordinates of the maximum of the 
normal moveout hyperbola are given by x m = —2d sin 9 and t m — to, one could, by making use of the standard 
trigonometric result, sin 2 9 + cos 2 9 = 1, derive the condition, 




The form of this conic section is that of an ellipse. From this it is now possible to conclude that the peak of the normal 
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Fig. 4. — For an inclined reflector, imaging the depth vertically from the surface of the earth, gives a false position of the reflector. 
The actual point of reflection on the inclined plane is indicated by N, but a vertical downward projection from S will give an 
erroneous impression that the reflection point is at D. The line SN is normal to the reflector interface. Migration will give a more 
correct impression of the true reflector interface. 



moveout curve will trace an elliptical trajectory in the x—t plane, with a continuous change in the dip angle, 9. 

While this entire algebraic exercise looks very elegant in principle, it fails hopelessly in addressing a serious 
practical problem. The derivation of Eq. (0]i has been done with the help of a single source-receiver pair. It assumes 
that every source-receiver pair in a full array of sources and receivers will have a common point of reflection from 
the inclined subsurface. A priori this is not possible to achieve in practice. If anything, for a realistic arrangement of 
source-receiver arrays, each source-receiver pair will have its own particular reflection point. Therefore, unlike the 
simple representation in Fig. Q] there will be no common depth point vertically below the common midpoint of the 
source receiver arrays. Being unmindful of this crucial point will only serve to furnish an erroneous impression of the 
subsurface depth, of the kind that has been schematically shown in Fig. [4] This fact categorically precludes the simple 
approach of determining the value of d, using the arguments which have been presented following Eq. (0, at the end 
of SectionQ] 

The difficulties do not end here. In actual fact the situation is vastly more complicated. The dipping reflectors 
will have a varying gradient from one spatial position to another, as opposed to being at a constant incline with respect 
to the horizontal, shown in Fig. [3] Besides this, there will be geological faults, and the propagation velocity of the 
seismic waves will also have a variation with respect to the depth. All of these facts will combine to render the accurate 
imaging of the subsurface a quite formidable task. Nevertheless, all reflection seismic records must be corrected for 
non-vertical reflections. This evidently complicated process of correction is what is, for all practical pur poses, implied 
by the term m igration. It purports to shift all dipping reflections to their true subsurface positions dLowrielll997 ; 
YilmazlboOll) . 



Another important point has t o be stressed here. So far the entire discussion has accounted for dip and subsurface 



variations in the inline direction (Yilmaz 2001), i.e. the direction along which the source-receiver arrays have been 
lined up. This means that the analytical treatment presented so far has been confined to a 2D plane only — along 
the vertical depth (given by the z coordinate) and along the offset direction (given by the x coo rdinate). How ever, 



practically speaking, variations in the subsurface will also take place in the crossline direction (Yilmaz 2001J), i.e. 
along the y axis. So a comprehensive migration process must account for all variations along all the three spatial 
coordinates, i.e. it should in effect be a complete 3D process accounting for any change in the gradient of any 
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subsurface reflector. 



3. Migration: A mathematical theory 



By now the case for a complete 3D migration process has been fully argued for. To carry out 3D migration 
effectively, a necessary exercise would be to study the propagation of a compressional wavefield, P = P(x, y, z, t). 
This is described mathematically by the wave equation 

1 d 2 P 



V 2 P 



dt 2 



0, 



(6) 



with V 2 being the standard Laplacian operator dArfken & Webei 2001), and with v having a general dependence on 
all the three spatial coordinates. Expressed in terms of its individual components, this second-order partial differential 
equation in both the space and the time coordinates will read as 

d 2 d 2 d 2 1 



dx 2 dy 2 dz 2 



. d 2 \ 

-5? P = - 



(7) 



On the surface of the earth the upcoming seismic wavefield would be measured as P(x, y, 0, t). For fixed values 
of x and y, the signal would be a function of t only, giving a single trace. Extended over the whole x-y plane, this will 
build a proper mathematical understanding of the entire 31? volume of seismic traces that one is required to work with. 
Given this boundary condition on the surface of the earth, one would have to find the wavefield, P(x, y, z, 0), which 
is actually the reflectivity of the earth in its recesses (and effectively, it also defines an initial condition at t — 0). This 
whole process will require, as an intermediate step, the projecting of the surface wavefield to a depth z, corresp onding 
to t = 0. In mathematical terms, this is the most general understanding that one can have about migration dYilmaz 



2001). 



Solving Eq. (|7]i is no easy feat, because of its dependence on four independent variables. However, under the 
simple wor king assumption that v is a constant, it is possible to go a long way, by adopting the method of Fourier 
transforms dArfken & Weberll2001). By definition, th e Fourier transform of the wavefield P(x, y, z, t) over the coor- 
dinates x, y and t is given as ( Arfken & WebeifeoOl ) 



P(k x ,k y , z,uj) 



P(x, y, z, t) exp (ik x x + ik y y — iuot) dx dy dt. 



(8) 



(2tt) 3 / 2 

In terms of this transformed wavefield, Eq. (O can be recast as a second-order ordinary differential equation, which is 



d 2 P 
dz 2 



k ) P — 0, 



(9) 



with P now being understood to be P = P(k x , k y , z, ui). 

It is worth stressing the importance of the method that has led to Eq. The starting point was Eq. © — 
a second-order partial differential equation in four variables. The Fourier transform of x, y and t, into their corre- 
sponding reciprocal domains of k x , k y and to, has rendered this partial differential equation into a relatively simple 
second-order ordinary differential equation. For mathematical ease, this is the most important purpose that a Fourier 
transform can serve. Of course, working in the Fourier domain has other practical advantages, especially from the 
perspective of the processing geophysicist, like the elimination of spikes and noise from the data (Marietal ■1 fl999l) . 



Following the Fourier transform, it will now be necessary to solve Eq. (O to extrapolate the wavefield from its 
value on the surface of the earth (at z = 0) to the required depth (where z ^ 0). Working under the simplifying 
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assumption that v is constant jYilmajboOlb. it becomes q uite easy to solve Eq. (O, which will act ually have th e form 
of the simple harmonic oscillator ( Arfken & WeberlboOl ). With the imposition of the definition dYilmazl200ll) 



k 2 - — 



(a dispersion relation as it wer^j), one can proceed to derive an integral solution of Eq. © as 

P(k x ,k y , z,uj) = P(k x ,k y ,0,u) exp(-ik z z) . 



(10) 



(ID 



Actually, the most general form of Eq. (fTTT > is given by the superposition of an upcoming wave and a downgoing 
wave. i.e. exp(±ifc z z). However, only the solution with the negative sign need be retained here, since it corre sponds 
physically to upcoming waves, which, for the present purpose, are actually the relevant modes ( Yilmaal200ll) . This 
whole analysis shows how it should be possible to extrapolate the surface wavefield to a given depth below the surface 
of the earth. While it has been quite a simple task to arrive at this result under the assumption of constant velocity, one 
can also derive an asymptotic solution for the extrapolation exercise when the velocity is dependent on the depth, i.e. 
v = v(z). This treatment has been presented in Section[4] 

One way or the other, the one primary fact that stands out is that so far as migration is concerned, the result given 
by Eq. (fTTT i. i.e. extrapolating the surface data to the depth, represents a crucial step. Beyond this point there can be 
many variations in the subsequent execution of the migration process, depending on varying contingencies. In what 
follows later, the specific principles of some of these methods will be dwelt upon. 



4. Extrapolating the wavefield in the case of depth-dependent velocity: The WKB approximation 



It has already been seen how easily Eq. (|9]l can be integrated when v is constant. Supposing now that v has a 
dependence on the z coordinate, i.e. v = v(z), it will still be possible to derive an asymptotic solut ion of Eq. ||9), under 
appropriate conditions , with the help of the WKB (Wentzel-Kramers-Brillouin) approximation ( Ma thews & Walker 
ll970l:lKing et al.ll2003b . 



Since v has a dependence on z, so should k z have a dependence on z, by virtue of Eq. ( [Tol l. Therefore, one can 



write 



d 2 P 

dz 2 



+ k 2 (z)P = Q, 



(12) 



for which a trial solution in the most general form can be set down as 

P = A{z)e la{z \ (13) 
with A being an amplitude function, and a being a phase function. Substituting this trial solution in Eq. (fT2l will give 



1 d 2 A 

A dz 2 



.da d 

dz dz 



ln ' d~z A 



0. 



(14) 



The real and imaginary components above will have to be set to zero separately. In doing this, a solution for the 
imaginary part can be extracted, and it will read as 



A=C 



da 



dz 



-1/2 



(15) 



1 In the case of the so-called Exploding Reflectors model, the velocity in Eq. (To) is transformed according to the rale, v ► v/% I Yilm azl200ll) . 
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where C is a constant of integration. The foregoing result indicates that knowing either A or a independently will 
suffice to give the full solution implied by Eq. (fT~3T > . 

So far the mathematical treatment has been exact. In deriving a solution from the real component of Eq. (TBI , it 
will be necessary now to apply the WKB approximation, requiring that the amplitude function A(z) in Eq. (fT3l l, is a 
much more slowly varying function of z than the phase, a(z). In that event one can neglect the second derivative of A 
with respect to z in Eq. ( TBI . This will give the phase solution 

(16) 



a(z) ~ — / k z (z)dz, 

(taking the more practically relevant negative solution, as usual) which, used along with Eq. <fT3T >. will give an aysmp 
totic solution that will read as 

C 



P 



: exp 



i / k z (z)dz 



(17) 



How can the use of the WKB approximation be justified self-consistently? On using the solutions given by 
Eqs. (Q3) and (fTol l. it can be seen that the second derivative of A is given as 



d 2 A 
ck 2 " 



.4 



u> 

k~v 



k,v 



-(Inv) 



d 2 

d^ (lnu) 



(18) 



On the right hand side of Eq. ([T8l one can see derivatives of the logarithm of v. Now it is known that the logarithm 
of any function varies much more slowly than that function itself. The same should hold for In v in comparison with 
v itself. Therefore, it becomes eminently possible to argue that with the second derivative of A depending on the 
derivatives of In v, the requirement of A being a slowly varying function of z is satisfactorily fulfilled, and one may 
safely neglect the term containing the second derivative of A in Eq. (TT~4b - Indeed, this will be all the more appropriate 
for greater depths, where, with v being an increasing function of z, the logarithmic terms will strongly temper the 
variation of A arising due to the velocity variations. So the asymptotic solution given by Eq. (T% will have a firmer 
validity at greater depths. This entire line of reasoning upholds the invoking of the WKB approximation. 



5. The ZD phase-shift migration in the frequency-wavenumber domain 



It has been discussed in Section[3]that the very definition of migration necessitates projecting seismic data from 
the boundary condition z — for all t, to the initial condition t = for all z. To this end, it shall be expedient to set 
down first the inverse Fourier transform corresponding to Eq. (|8). This will be 

1 



P{x,y,z,t) 



(2tt) 3 / 2 



P(k x , k y , z, to) exp {—ik x x — ik y y + iujt) dk x dk y du. 



(19) 



It is also known by now how it should be possible to relate P(k x , k y , z, ui) to the surface seismic data with the help of 
Eq. (fTTTl . Using this condition in Eq. (fT9l and setting t = Q, will give 

1 



P(x,y,z,0) 



(2tt) 3 / 2 



P{k x , k y , 0, w) exp (—ik x x — ik y y — ik z z) dk x dk y dw, 



(20) 



which is the equation for the phase-shift method, suited for 3D migration (lYilmaauOOll) . The involvement of u> can 
be eliminated from thi s equation, b y making use of Eq. ( fTUl ), under the constraint that v is constant, and k x and k y are 
to be kept unchanged ( Yilmaz 20011) . This will first deliver the relation 

vk z 



dw = 



1.2 I UI I U2 
rh x T hjy ^ tv z 



--dk z 



(21) 
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leading ultimately to 

P(x,y,z,0) 



vk. 



( 2 ^) 3/2 J J J Jki + kj+ ki 



P [k x , ky , 0, V\ k x -\- ky + k z 



x exp {—ik x x — iky-y — ik z z) dk x dk y dk z 
which is the equation for constant velocity 3-D Stolt migration jYilmazlboOlT) . 



(22) 



The derivation of the result in Eq. (1221 is certainly consistent with what has been forwarded as a mathematical 
definition of migration in Section [3] However, short of resorting to an involved numerical treatment, the full 3D 
problem does not lend itself very easily to deriving any pedagogical insight out of it. Hence it should be instructive to 
consider a very special case of this migration method, for a wavefield that has a spatial dependence on the z coordinate 
only, and so will actually be a vertically propagating planar wavefield. 

For such a wavefield, with no dependence on the x and y coordinates, and propagating in a planar front along 
the z axis only, one can set down the wavefield as P = P(z, t). In terms of the Fourier inverse transform it can be 
expressed as 



P(z,t) 



1 



P(z, lo) exp {—iujt) doj. 



(2tt) 1 /2 

Now it is already known that Eq. (fTTT i will give an extrapolation relation going as 

P{z, uj) = P(0, oS) exp (—ik z z) , 



(23) 



(24) 



for the special case that is being studied here. From Eq. (TT~0b one also gets k z = oj/v, since there is no involvement of 
the x and y coordinates. Use of all these conditions in Eq. ( f23l , will result in 



P(z,t) = 



(2tt) 1 /2 



P(z, to) exp 



+ t 



which, for t = 0, will give 



P(z,0) 



(2tt) 1 /2 



P(0, to) exp 



-ioj 



(- 

V V 



du), 



duo = P(0,z/v). 



(25) 



(26) 



But it is also known that t and ui are conjugate variables for the Fourier transform. Their mathematical connection is 
given by 



P(0,t) 



1 



(2tt) 1 /2 

which can then be compared with the result given in Eq 



P(0,uj)exp(-iujt)duj, (27) 
26]). This will immediately allow for setting down 

z = vt, (28) 



which is a result that very unambiguously connects the depth of a reflecting stratum with the time of the reception of a 
signal reflected from that stratum (under the simple condition that the velocity of the wave propagation is unchanging). 
An important point to note here is that the simple appearance of Eq. d28l has been possible because the plane wave 
allows only the vertically propagating mode in Eq. ( [Tol l, and so by default this mode has been decoupled from all the 
other modes. This decoupling should be impossible to achieve for higher spatial dimensions, as a look at Eq. (TlOb 
ought to reveal. 
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6. The stationary phase approximation 

It shall be useful at this point to derive the traveltime equation for the 3-D zero-offset case, using the method of the 
stationary phase approximation. From Eq. ( fT9l ) one can, making use of the extrapolation condition given by Eq. dTTb , 
write a compact relation that will read as 

P(x,y,z,t) = , 2 \ 3/2 J J J P(k x ,ky,0,w)exp(i$z)dk x dkydu, (29) 

in which 

k x x k v y tut 

$ = -k z - — v — + — . (30) 

z z z 



It is easy to see that the integrand in Eq. d29j l has an amplitude part and a phase part. If the phase is to vary much 
more rapidly than the amplitude (a requirement that is consistent with the WKB analysis presented in SectionHJ), then 
the consequent rapid oscillations over most of the range of the integrati on, will result in an average value of nearly 



zero. The exception to this principle will occur when <£> is stationary (Ma thews & Walked 1 1970t lArfken & Webei 



2001). This (the stationary phase approximation) will be quantified by i5$ = 0, which is, in fact, a condition for 
the turning points of $. The major contribution to the integral will, therefore, come from the turning points of the 
phase function $. This condition will help in identifying the physical coordinates from where one will obtain the most 
significant contribution to the signal received on the surface of the earth. 

Going by Eq. ( fTOt , it is already known that k z has a dependence on k x , k y and u. Taking this in conjunction with 
Eq. d30l ), and noting that the relevant variables for the integration in Eq. d29l are k x , k y and ui, the required extremum 
condition for $ can consequently be stated as 

s *={wJ §k * + {wJ Sk « + {^) 5 " = °- (31) 

The result above can only mean that all the three partial derivatives in Eq. d3"T1 > are to be individually set to zero. This 
will give three distinct mathematical conditions, from which one will be able to derive the traveltime equation for the 
3D zero-offset case as 

x 2 + y 2 + z 2 = {vt) 2 . (32) 

This defines a sphere of radius vt in the x-y-z coordinate space, and a hyperboloid in the x-y-t space for a constant 
value of z. Considering only the inline direction and the depth (i.e. when y = 0), Eq. (l32l will define a semi-circle 
below the surface of the earth in the x—z plane, and the usual hyperbola in the x—t plane. 

However, a much more interesting result can be seen for the simplest possible case of the planar wavefront, i.e 
where the spatial dependence is only on the z coordinate. In that case the result in Eq. (l32l will converge simply to the 
linear solution z = vt, which is exactly what Eq. d28l > also furnishes. So effectively this simple special case establishes 
an intriguing mathematical equivalence between the stationary phase approximation method and the constant velocity 
Stolt migration algorithm. 



7. Some general comments 

Migration strategies can be implemented over a wide range — from 2D poststack time migration to 3-D prestack 
depth migration. The geological features of the subsurface will determine the choice of a particular algorithm, or even 
a combination of various algorithms. In practice, time migration enjoys wide popularity because it is least sensitive 
to errors in velocity, and delivers acceptable results to make a reliable interpretation. However, time migration is 
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trustworthy only as long as the lateral velocity variations are at most moderate. Otherwise depth migration is more 
appropriate. 

Standard migra tion algori thms are based on the wave equation. One can classify migration processes under three 



broad categories (lYilmazE oOl) 



• Those based on the integral solution of the wave equation. 

• Those based on the finite-differencing of the partial differential equation implied by the wave equation. 

• Those based on the implementation of the frequency-wavenumber algorithm. 

It is to be borne in mind that none of the migration algorithms that are implemented in the industry today gives 
a complete prescription for handling all dips and velocity variations, while being cost-effective at the same time. Mi- 
gration algorithms based on the integral solution of the wave equation (Kirchhoff migration) can become cumbersome 
when they encounter lateral velocity variations. Finite-difference algorithms can only be specified for different degrees 
of dip approximations. And frequency-wavenumber algorithms are effective against lateral velocity variations only to 
a limited extent. 
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